function lines=div_lines(img,size,sigma,is_show)
    if ~exist('lines.mat','file')
        % img: input graphs, should be m x n matrix instead of tensor!!
        %[gx,gy,gxx,gyy,gxy]=prepare(img,size,sigma);
        %% Divide each box using DBSCAN Cluster method
        img=255-img';
        bool=0*img;
        bool(find(img>30))=255;
        [row,col]=find(bool>0);
        boxes=[[row],[col]];
        idx = dbscan(boxes,4,2);
        %% For each box,using GMM method to divide lines
        % those who cannot be divided will be blocked in output tuple. 
        num_ind=max(idx);
        lines.box={};
        lines.idx={};
        lines.valid={};
        lines.length=num_ind;
        for i=1:num_ind
            box=boxes(find(idx==i),:);
            if length(box)>50000 || length(box)<100
                % filter those who owns too few and too many points.
    %             boxes(find(idx==i),:)=[];
    %             idx(find(idx==i),:)=[];
                lines.valid{i}=0;
                continue
            end
            %gscatter(box(:,1),box(:,2),idx(find(idx==i),:));
            try
                GMModel = fitgmdist(box,6,'Options',statset('TolFun',1e-5));
                idx2=cluster(GMModel,box);
                %gscatter(box(:,1),box(:,2),idx2);
                lines.box{i}=box;
                lines.idx{i}=idx2;
                lines.valid{i}=1;
            catch
                lines.valid{i}=0;
                continue
            end
        end
        save('lines.mat','lines');
    else
        load('lines.mat','lines');           
    end
    if is_show==1
        for i=1:lines.length
            if ~lines.valid{i}
                continue 
            else
                gscatter(lines.box{i}(:,1),lines.box{i}(:,2),lines.idx{i});
                hold on;
                text(mean(lines.box{i}(:,1)),mean(lines.box{i}(:,2)),int2str(i));
                hold on;
            end
        end 
    end
end